Référence pour les séries chronologique :

https://otexts.com/fpp2/

Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. OTexts.com/fpp2. Accessed on 2 may 2022

Le jeu de données

Source : enquête fréquentation touristiques (https://www.ispf.pf/bases/Tourisme/EFT.aspx)

library(data.table)
library(ggplot2)
library(hts)
library(plotly)

data = fread("SerieTourisme.csv")

# la structure
str(data)
## Classes 'data.table' and 'data.frame':   8725 obs. of  4 variables:
##  $ date     : IDate, format: "2006-01-01" "2006-02-01" ...
##  $ Region   : chr  "" "" "" "" ...
##  $ PaysDiff : chr  "" "" "" "" ...
##  $ Touristes: int  108 56 47 46 54 24 21 41 21 34 ...
##  - attr(*, ".internal.selfref")=<externalptr>
# quelques statistiques
data[,Region:=as.factor(fifelse(Region=="",NA_character_,Region))]
data[,PaysDiff:=as.factor(fifelse(PaysDiff=="",NA_character_,PaysDiff))]
summary(data)
##       date                             Region             PaysDiff   
##  Min.   :2006-01-01   Europe (hors France):2993   Allemagne   : 192  
##  1st Qu.:2009-10-01   Asie                :2113   Autre Europe: 192  
##  Median :2013-09-01   Pacifique           :1336   Belgique    : 192  
##  Mean   :2013-09-26   Amérique du Sud     : 742   Canada      : 192  
##  3rd Qu.:2017-08-01   Amérique du Nord    : 573   France      : 192  
##  Max.   :2022-03-01   (Other)             : 956   (Other)     :7753  
##                       NA's                :  12   NA's        :  12  
##    Touristes     
##  Min.   :   1.0  
##  1st Qu.:  14.0  
##  Median :  43.0  
##  Mean   : 326.1  
##  3rd Qu.: 194.0  
##  Max.   :8748.0  
## 
# la série sous format TS 
Touristes    <- data[,lapply(.SD,sum,na.rm=T),.SDcols="Touristes", by="date"]
Touristes.ts <- ts(Touristes$Touristes, 
                   frequency = 12, 
                   start = Touristes[date==min(date), c(year(date),month(date))],
                   end = Touristes[date==max(date), c(year(date),month(date))])
print(Touristes.ts)
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2006 13714 15176 16256 18061 18051 18138 22917 20035 21144 21120 18401 18536
## 2007 15354 15859 18107 17290 17818 18783 21034 20897 18789 20891 16229 17190
## 2008 13102 14765 16829 15962 17485 16551 19111 18601 17989 18121 14040 13940
## 2009  9999 10372 12415 11230 13236 13824 16853 15808 14888 15972 12892 12958
## 2010  9016  9730 10547 10271 11525 12119 17790 15087 15160 16092 12784 13798
## 2011 11371 11038 12304 12458 12838 14424 16858 15372 14402 14519 13086 14106
## 2012 10238 11523 13075 13147 13879 14940 16979 16002 15944 15519 12470 15262
## 2013 11174 11177 13897 12011 13534 15120 17289 14655 14175 14576 12953 13832
## 2014 12422 12410 15410 15737 14853 14650 17656 14603 15500 17546 14646 15169
## 2015 12343 12949 14472 13956 14832 16223 18060 16463 16927 17561 15681 14364
## 2016 12340 14444 15168 17042 15738 16758 19517 17719 17117 17499 14295 14858
## 2017 11910 13564 16281 15523 16782 17599 21448 18563 18539 18271 15449 15030
## 2018 11457 15747 17452 15955 16559 19372 24168 20110 19809 20661 17241 17737
## 2019 15007 16752 18684 19240 18749 21487 25361 21864 20305 21170 19185 18838
## 2020 13948 15497  7491  4605  7834  7680  8976  4486  6500  3924   524   347
## 2021  3368  8552 14331  9481  6739 12346 10320 12321  6262  8982 14856   293
## 2022 13714 15176 16256
derAnnee <- floor(max(time(Touristes.ts)))
derMois  <- round((max(time(Touristes.ts)) %% 1)*12 + 1)

ggplotly(
  autoplot(Touristes.ts, ylab="", xlab="", main="Evolution du nombre de touristes") +
    geom_point(data=window(Touristes.ts,c(derAnnee,derMois),c(derAnnee,derMois)), size=2,color="red",show.legend=F))

Etude de la saisonnalité

Brute

ggplotly(
  ggseasonplot(Touristes.ts, ylab="", xlab="") + 
    ggtitle("Evolution du nombre de touristes par mois") +
    guides(colour=guide_legend(title="Années",ncol=2)))
ggplotly(
  ggmonthplot(Touristes.ts, ylab="", xlab="") +
    ggtitle("Moyenne des variations annuelles par mois du nombre de touristes"))

Sur fenêtre de temps

FirstAnnee    <- floor(min(time(Touristes.ts)))
FirstMois     <- round((min(time(Touristes.ts)) %% 1)*12 + 1)
LastAnnee     <- 2019
LastMois      <- 12
Touristes.ts2 <- window(Touristes.ts,c(FirstAnnee,FirstMois),c(LastAnnee,LastMois))

ggplotly(
  ggseasonplot(Touristes.ts2, ylab="", xlab="") + 
  ggtitle("Evolution du nombre de touristes par mois") +
  guides(colour=guide_legend(title="Années",ncol=2)))
ggplotly(
  ggmonthplot(Touristes.ts2, ylab="", xlab="") +
  ggtitle("Moyenne des variations annuelles par mois du nombre de touristes"))

Décompostion de la série

L’hypothèse retenue sur la saisonnalité est qu’elle ne varie pas au cours du temps. Le résidus est le résultat de la série d’origine à laquelle on soustrait la saisonnalité et la tendance estimées.

Décompostion additive, multiplicative (classique)

library(dplyr)
ggplotly(
  Touristes.ts2 %>% decompose(type="additive") %>%
    autoplot() + xlab("Année") +
    ggtitle("Décomposition saisonnière classique")
)

Seasonal and trend decomposition using Loess STL (modèle par régression locale)

Les fonctions à retenir :
* stl()
* seasonal()
* trendcycle()
* remainder()
* seasadj ()

Touristes.stl <- stl(Touristes.ts2, s.window="periodic", robust=TRUE) 

remain.ts     <- remainder(Touristes.stl)
remain.temps  <- time(remain.ts)[which(abs(remain.ts)>quantile(abs(remain.ts),probs=0.95))]
remain.valeur <- remain.ts[which(abs(remain.ts)>quantile(abs(remain.ts),probs=0.95))]
remain.dt     <- data.table(datetime=remain.temps,y=remain.valeur,parts=factor("remainder"))

ggplotly(
  autoplot(Touristes.stl) + ggtitle("Décomposition de la série (data = seasonal + trend + remainder)") +
    geom_point(data=remain.dt,aes(x=datetime, y=y), size=2,color="red",show.legend=F) +
    xlab("") + ylab("Touristes")
)

Série CVS (Corrigée des Variations Saisonnières)

La fonction seasadj peut prendre différents objets TS, comme la décomposition classique

ggplotly(
  autoplot(cbind(brute=Touristes.ts2,CVS=seasadj(Touristes.stl)), ylab="Touristes") + 
  ggtitle("Série corrigée des variations saisonnières") + 
  scale_colour_manual(values = c("dark gray","blue")) +
  xlab("")+
  theme_bw() + 
  theme(legend.position = "bottom")
)

Analyse multivariée par décompostion hiérarchique de la série

TS multivariées

# les séries sous format TS 
data[,Region2:=Region]
data[Region %in% c("France","Europe (hors France)"),Region2:="Europe"]
PaysImp <- data[,.(Touristes=sum(Touristes)),.(PaysDiff,Region2)][Touristes>=50000, .(PaysDiff,Region2)]
data[, c("Region3","PaysDiff2"):=.("Autre","Autre")]
data[PaysImp, c("Region3","PaysDiff2"):=.(i.Region2,i.PaysDiff), on=.(Region2,PaysDiff)]
data[,c("Region3","PaysDiff2"):=.(sprintf("%10s",abbreviate(Region3,10)),
                                  sprintf("%7s",abbreviate(PaysDiff2,7)))]
print(data[,.(Touristes=sum(Touristes)),.(Region3,PaysDiff2)])
##        Region3 PaysDiff2 Touristes
##  1:      Autre     Autre    370260
##  2: AmériqdNrd    Canada    103476
##  3: AmériqdNrd       USA    918050
##  4:       Asie     Japon    201141
##  5:     Europe   Allemgn     60058
##  6:     Europe   Espagne     51547
##  7:     Europe    Italie    142705
##  8:     Europe   Roym-Un     61146
##  9:     Europe    France    640486
## 10:  Pacifique   Austral    131777
## 11:  Pacifique   Nvll-Cl     62953
## 12:  Pacifique   Nvll-Zl    101660
TouristesParPays    <- dcast(data, date ~ Region3 + PaysDiff2, fun.aggregate = sum, value.var = "Touristes",sep="")
TouristesParPays.ts <- ts(TouristesParPays[,-1], 
                          frequency = 12,
                          start = TouristesParPays[date==min(date), c(year(date),month(date))],
                          end = TouristesParPays[date==max(date), c(year(date),month(date))])

ggplotly(autoplot(TouristesParPays.ts, ylab="", xlab="",colour = T,facets = T,
                  main="Evolution du nombre de touristes par pays",
                  space="free_y") + theme(strip.text.y = element_blank()))

HTS (série hiérarchique)

TouristesParPays.hts <- hts(TouristesParPays.ts, characters = c(10,7))
str(TouristesParPays.hts)
## List of 3
##  $ bts   : Time-Series [1:195, 1:12] from 2006 to 2022: 1599 1724 1758 1528 1745 1630 1998 1796 2117 2030 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:12] "      Asie  Japon" "     Autre  Autre" "    Europe France" "    Europe Italie" ...
##  $ nodes :List of 2
##   ..$ Level 1: int 5
##   ..$ Level 2: 'table' int [1:5(1d)] 1 1 5 3 2
##   .. ..- attr(*, "dimnames")=List of 1
##   .. .. ..$ : chr [1:5] "      Asie" "     Autre" "    Europe" " Pacifique" ...
##  $ labels:List of 3
##   ..$ Level 0: chr "Total"
##   ..$ Level 1: chr [1:5] "      Asie" "     Autre" "    Europe" " Pacifique" ...
##   ..$ Level 2: chr [1:12] "      Asie  Japon" "     Autre  Autre" "    Europe France" "    Europe Italie" ...
##  - attr(*, "class")= chr [1:2] "hts" "gts"
ggplotly(autoplot(aggts(TouristesParPays.hts,levels = 0), ylab="", xlab=""))
ggplotly(autoplot(aggts(TouristesParPays.hts,levels = 1), ylab="", xlab="",colour = T,facets = T) + theme(legend.position="none"))
TouristesParPays1.stl.list   <- lapply(aggts(TouristesParPays.hts,levels = 1), FUN = function(x) stl(x, s.window = "periodic", robust = T))
TouristesParPays1.trend.list <- lapply(TouristesParPays1.stl.list, FUN = function(x) trendcycle(x))

ggplotly(autoplot(
  ts(as.data.table(TouristesParPays1.trend.list),
     frequency = 12,
     start = TouristesParPays[date==min(date), c(year(date),month(date))],
     end = TouristesParPays[date==max(date), c(year(date),month(date))]),
  main="Tendance du nombre de touristes par région", ylab="", xlab=""))

Projection du nombre de touristes

Projection selon un modèle ARIMA

Projections sur les 12 mois à venir issues d’un modèle ARIMA. Le modèle a été élaboré suite à l’analyse des graphes d’autocorrélation et d’autocorrélation partielle. Le modèle a été testé sur plusieurs échantillons tests et mis en compétition avec un modèle de sélection automatique des paramètres.

# Autocorrélation et Autocorrélation partiel
Touristes.ts2 %>% ggtsdisplay(main="")

# Autocorrélation et Autocorrélation partiel sur série désaisonnalisée
Touristes.ts2 %>% diff(lag=12) %>% ggtsdisplay(main="Série désaisonnalisée")

# Autocorrélation et Autocorrélation partiel sur série désaisonnalisée et sans tendance
Touristes.ts2 %>% diff(lag=12) %>% diff() %>% ggtsdisplay(main="Série désaisonnalisée et sans tendance")

# Modèle arima
fitarima <- arima(Touristes.ts2, order=c(0,1,4), seasonal = list(order=c(0,1,1), period=12)) 
checkresiduals(fitarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,4)(0,1,1)[12]
## Q* = 19.803, df = 19, p-value = 0.4065
## 
## Model df: 5.   Total lags used: 24
# AIC pour comparer les modèles
AIC(fitarima)
## [1] 2586.577
# Projection
fc       <- forecast(fitarima, h=12 , level=c(95))

autoplot(fc) +
  ggtitle("Projection du nombre de touristes selon un modèle ARIMA",
          subtitle = "Paramètres du modèle p=0, d=1, q=4, P=0, D=1, Q=1, s=12") +
  guides(fill=guide_legend(title = "Intervalle \nde confiance" )) +
  ylab("") + 
  xlab("")

Projection hiérarchique

Projections sur les 12 mois à venir issues d’un modèle hierarchique à 2 niveaux (pays et régions). Pour chaque série de la hierarchie les projections de base sont obtenues à l’aide d’un modèle automatique ARIMA. Les projections de base sont révisées de telle sorte que les aggrégations par niveau soient conséquentes. Par exemple, l’aggrégat des projections des pays d’une région correspond aux projections de la région en question.

library(parallel)
TouristesParPays2.hts <- hts(window(TouristesParPays.ts,c(FirstAnnee,FirstMois),c(LastAnnee,LastMois)) , characters = c(10,7))
fc.hts                <- forecast(TouristesParPays2.hts, h=12, fmethod="arima", parallel = T, num.cores=detectCores())

ggplotly(
  autoplot(allts(fc.hts, forecasts = T), ylab="", xlab="")
  )